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ABSTRACT 

V 

A  fast  algorithm  for  reconstructing  images  governed  by 
a  2-D  Helmholtz  equation  is  presented.  The  computa¬ 
tional  complexity  of  the  algorithm  is 
O(iWlogM)  or  0(iVM2)  depending  on  boundary  condi¬ 
tions,  where  N  and  M  are  the  number  of  spatial  grid 
points  in  the  x  and  y  directions  respectively.  This  prob¬ 
lem  arises  when  smoothing  a  large  number  of  images 
governed  by  the  2-D  wave  equation,  because  a  Fourier 
transform  in  time  gives  a  new  set  of  images  governed  by 
the  Helmholtz  equation.  When  the  images  come  from  a 
scattering  process,  we  show  that  a  linear  least-squares 
Born  inversion  of  the  wave  field  amplitudes  can  be  per¬ 
formed  during  the  smoothing  procedure  without  chang¬ 
ing  the  computational  complexity.  We  also  show  that 
the  smoothing  algorithm  is  well-posed,  and  that,the  sam¬ 
ple  functions  of  the  smoothed  estimate  possess  smooth¬ 
ness  properties  consistent  with  the  Helmholtz  equation.  _ 


1.  Introduction 

In  this  paper  we  derive  a  fast,  recursive,  linear 
least-squares  smoothing  algorithm  for  the  2-D 
Helmholtz  equation.  Our  algorithm  can  be  used,  for 


example,  to  smooth  a  large  number  of  images  governed 
by  the  2-D  wave  equation  arising  in  acoustical  hologra¬ 
phy  [8]  or  in  oceanic  surveillance.  If  we  assume  that  the 
vibrating  system  is  in  steady  state,  and  that  the  inputs 
and  observation  noise  are  temporally  stationary,  then  a 


*  This  work  was  supported  by  the  Office  of  Naval  Research  under  Contract  N00014-85-K-0255 


C-Xv IV 


Fourier  expansion  with  respect  to  time  gives  a  new  set 
of  images,  indexed  by  the  frequency  /,  that  are 
governed  by  the  Helmholtz  equation  whose 
wavenumber  k  varies  with  /.  Each  transformed  image 
can  be  smoothed  separately,  and  estimates  of  the  origi¬ 
nal  time-domain  images  can  then  be  obtained  by  an 
inverse  Fourier  transform. 

In  order  to  smooth  images  governed  by  the 
Helmholtz  equation,  we  reformulate  the  equation  as  a 
well-posed,  distributed-parameter,  acausal  linear  system, 
and  use  the  recent  extension  of  Adams,  Willsky  and 
Levy  [I]  of  the  method  of  complementary  models  [ 9 j  to 
write  the  Hamiltonian  system  for  the  smoothed  estimate. 
Transforming  in  one  direction  produces  a  set  of  indexed, 
well-posed,  finite-dimensional,  acausal  linear  systems 
which  we  solve  recursively  using  a  diagonalizing  change 
of  variables.  The  complexity  of  our  algorithm  for  each 
wavenumber  k  is  0(NM\ogM)  or  0(NM2)  depending  on 
the  boundary  conditions,  where  N  is  the  number  of  grid 
points  in  the  x  direction  and  M  is  the  number  in  the 
(orthogonal)  y  direction. 

Another  approach  to  the  smoothing  of  images 
governed  by  the  Helmholtz  equation  inyolves  the  use  of 
the  Karhunen-Loeve  expansion  of  the  wave  field 
[2 j , [3] , [4] ,  resulting  in  an  algorithm  with  complexity 
0(MN\ogMN).  However,  the  boundary  conditions  are 
required  to  be  conservative  and  have  no  random  inputs. 
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In  contrast,  we  can  handle  x-boundary  conditions  that 
are  conservative  or  dissipative  and  that  include  random 
inputs.  We  note  that  work  using  a  similar  approach  is 
discussed  by  Yoshida  and  Ogura  in  [10].  In  their  work, 
the  dynamics  are  discrete  and  the  underlying  random 
field  is  homogeneous  (spatially  stationary),  whereas  in 
this  ,  the  dynamics  are  continuous  and  the 

random  field  is  not  required  to  be  homogeneous.  Furth¬ 
ermore,  an  important  step  in  the  derivation  of  the  esti¬ 
mator  in  [10]  is  the  replacement  of  a  non-Markovian 
random  process  with  a  Markovian  random  process  hav¬ 
ing  the  same  mean  and  covariance.  In  our  approach, 
this  realization  step  is  not  needed. 

2.  Problem  Statement 

Consider  the  scalar  Helmholtz  equation  on  the  rec¬ 
tangle  n  =  [0,L d  x  [0,Z,2]  : 

ux  z+Uyy+k-u  =  e{x,y)  (2.1) 

u  =  u(x,y)  ,  (x,y)  €  fi 
with  boundary  conditions 

u(0,y)  =  v{(y)  ,  u(Ll,y)  =  v 2(y)  (2.2a) 

«(x,0)  =  0  ,  u(x,L2)  =  0  •  (2.2b) 

Here 

k2  =  k£+jri  ,  k0  =  ,  n>  0 

c 

f  is  the  input  field,  iq  and  v2  are  boundary  inputs  on  the 
x-axis,  u  is  the  wavefield  amplitude,  /  is  the  temporal 


frequency  of  interest,  c  is  the  phase  velocity  of  the 
medium,  and  r?  is  the  damping  term.  The  observations 


z(x, y)  =  u(x,y)+w(x,y)  ,  ( x,y )  6  Q 
where  w(x,y )  is  the  observation  noise  field. 


(2.3) 


We  shall  make  the  following  statistical  assumptions: 
(1)  the  driving  field  e  and  observation  noise  field  w  are 
zero  mean  and  white  with  constant  intensities  q  and  r, 
respectively,  and  are  uncorrelated  with  each  other,  (2)  if 
v{y)  =  [v i(j/)  t>2(y)]' then 

Ev(y)  =  0 

Ev(y)v(s)'  =  Ylv6(y-s) 

where  n„  is  invertible  and  v  is  uncorrelated  with  e  and 

w. 

The  estimation  problem  of  interest  here  is  to  deter¬ 
mine  the  linear  least-squares  estimate  u(x,y)  of  u(x,y )  , 
(j ,y)  e  H,  given  the  observations  (2.3)  over  the  entire 
rectangle  n  . 

To  see  how  the  Helmholtz  equation  may  arise  in 
practice,  consider  the  2-D  wave  equation  on  the  rectan¬ 
gle  n  =  [0 ,L  i)  X  [O.Loj: 

c'2{uzz+uyy)+lut  =  d{x-y<t) 

u  =  u(x,y,t)  ,  ( x,y,t )  6  H  X  [T0,r,] 
with  boundary  conditions 


1  •*.  •*,  •*»  '  •  **•> 


M(0 ,y,t)  =  Vi(!/,t)  ,  u(Lvy,t)  =  v2(y,t ) 
u(x,0,t)  =  0  ,  u(x,Lo,t)  =  0 

u(x,y,T0 )  =  ut(x,y,T0)  =  0 
and  observations 

z{x,y,t)  =  u(x,y,t)+w(x,y,t) 

Let  T0— *  -  oo  and  assume  that  the  observation  interval  is 
[0, Tx]  where  Tx  is  very  large.  Then  a  Fourier  series 
expansion  of  the  observations  will  give  a  new  set  of 
images  z(x,y,f)  that  are  governed  by  (2.1)-(2.3),  where 
the  dependence  on  /  has  been  suppressed,  and  where 
e  =  -  d/c2  and  rj  =  271-/7.  We  assume  that  the  input  field 
d(x,y,t),  boundary  inputs  vx(y,t ),  v2(y,t),  and  observation 
noise  w(x,y,t)  are  wide-sense  stationary  in  time.  The  sta¬ 
tionary  assumption  implies  that  for  /it^/2  »  z{x,y,f i) 
and  u(x,y,fl )  are  uncorrelated  with  z[x,y,f2),  since 
T^  —  oo  .  We  can  therefore  solve  an  uncoupled 

set  of  smoothing  problems  for  the  2-D  Helmholtz  equa¬ 
tion  (indexed  by  /),  and  then  inverse  transform  to 
recover  the  time-domain  estimates. 

3.  State  Space  Formulation  and  Characterization  of 
the  Estimate 

In  order  to  put  (2.1)-(2.3)  in  state-variable  form, 
define  an  operator  T  with  domain  D(T )  as  follows: 


Tf  =-(k2f(x,y)+fyJx,y)),  f  6  D(T) 


where 


¥ 


D\T)=\f  €  Lo(n  ):/,/y  abs.  cont. ,  /yy  6  L2(fi  ), 

/  ( x  ,0)  =  f{x,L2 )  =  0} 

Also  define  the  state  vector  m(x,y)  as 

m(  x,y)  =  [m,(i,y),m2(i,y)]'  =  [u(x  ,y),ux{x  ,y)]' 
We  can  now  rewrite  ( 2.1)-( 2.3)  as 

Q 

-—m(x,y)  =  Am(x,y)+B(.(x,y) 


(3.1a) 


where 


ml  e  D(T) 


v(y )  =  V0m(0,y)  +  VL  m(Lvy) 


z(x,y)  =  Cm(x,y)+w(x,y ) 


-TO  ’  B  ~~  l]  ’  C  ~  [*  °]  ’  U  ~  v[ 


(3.1b) 

(3.1c) 


y  =  vr  — 

0  00  ’  VL  —  10 


Equations  (3.1)  are  in  the  form  of  a  distributed  parame¬ 
ter,  acausal  linear  system.  Finite-dimensional  acausal 
linear  systems  are  discussed  in  [6]-[7].  We  show  in 
Appendix  A  that  (3.1)  are  well  posed,  in  the  sense  that 
if  e  =  v  =  0,  then  m  =  0. 

Using  results  in  [l],  we  can  express  the  linear  leasts 
squares  estimate  m  of  m,  given  the  observations  (2.3), 
as  the  solution  of  the  following  Hamiltonian  system: 


x,y)  _  A  m{x,y) 


x,y)  r'xC'C  -A'  \\\(x,y)  - C'r-lz(x,y ) 


(3.2a) 


- 


rii  £  D (A)  ,  \  (E  D (A  ) 


0  = 


where 


v*n;lv 


m(0,y) 

4- 

o' 

\ko,y) 

m(Lvy) 

I 

0  I 

>> 

, 

(3.2b) 


^  =  VL\ 

the  asterisk  denotes  adjoint,  and  D(A)  ,  D(A*)  denote 
the  domains  of  A  and  A*.  These  domains  are 

D(A)  =  {[/ 2] ,;/i  €  D{T)  ,  f2  6L2(fl)} 


D(A*)  =  {[fJzY-.fx  6  L2(Q  )  ,  /2  6  D(T)} 

Moreover,  the  input  estimate  satisfies 

t  (x,y)  =  qB*\(x,y)  (3.2c) 

Eq.  (3.2c)  can  be  interpreted  as  a  linear  least-squares 
Born  inversion  when  the  observations  are  the 

scattered  wavefield  in  an  inverse  scattering  experiment. 

Instead  of  solving  (3.2)  directly  for  m  and  X,  we  will 
transform  (3.2)  with  respect  to  y  using  the  discrete  sine 
transform  S,  given  by 

h 

s 9  =  —  f  s\n{py)g{y)dy 
Li  0 


2nn  n  .  ,  ,  0 

p  =  -7 -  ,n  =  0,±1,±2, 


It  can  easily  be  verified  that  ST  =  (p2-  k2)S  and  thus 


SA  =  where 


A,  = 


0  1 

o  ,  o  _ 

p~- k~  0 


/.V.V.A 
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Transforming  (3.2a)  with  respect  to  y  gives 

o  I" „  ,  .  "I  A.  qBB'  %  i  r  „ 

J_rn(x.p)  p  m(x,p)  0  /ooflN 

dx  \(x,p)  r~1C'C  -A*  X( x ,p)  -C'r~1z(x,p)  ^  '  ’ 

where  m(x,p )  =  S m(x,y)  ,etc  .  Transforming  (3.2b)  with 


respect  to  y  gives 


n  —  V*T\~^V  ™(0iP)  4.  ^  0  ,X(0,p) 

U  V  1 A  fi  V  -  /  r  \  ^  n  T  \  /  T  \ 


(L  i,p)  T  0  /  X(L1)P)  • 


When  written  in  the  standard  form  for  an  acausal  linear 
system,  these  boundary  conditions  are 


0 Wn  ^(Tp)  ,  W  ™(Li,p) 

°  0  X(0,P)  +  ^  X(L  j,p) 


where 


(3.3b) 


^0  = 


Vo nv- 1 K0  -/  0 

^•n,-'F0  0  “dvVi=  n'n.-'n  / 


We  see  that  in  the  p  domain  the  estimate  Hamil¬ 
tonian  is  decoupled;  that  is,  (3.3)  are  simply  indexed  by 
p.  and  can  be  solved  individually.  In  the  next  section 
we  show  how  to  solve  (3.3)  recursively  for  each  p. 
Before  doing  so,  we  first  consider  whether  (3.3)  is  well- 
posed.  These  equations  can  be  shown  to  be  well-posed 
by  first  realizing  that  for  each  p,  (3.3)  is  the  estimate 
Hamiltonian  of  another  acausal  linear  system,  the  so- 
called  p-dynamics  of  the  Helmholtz  equation. 


—  m{X'p)  =  Apm(x  ,p)  +Bt(x  .p ) 


V'0m(0.p)+  VL  m(L  ,,p)  =  v(p) 


(3,1a) 


(3,1b) 


^ -v ■f ly t ^ -or  ..y  -yv '■  o  o ‘  -■ 


c  ( x ,  p )  =  Cm  (x  ,p)  +  w(x  ,p) 


(  3.1c) 


where  we  have  transformed  (3.1)  into  the  p-domain  by 
S.  The  in vertibility  of  V0  +  VL  eA,Ll  is  necessary  and 
sufficient  for  (3.4)  to  be  well-posed  [6] -[7].  It  is  easy  to 
verify  that  this  matrix  is  invertible  for  all  p.  To  prove 
then  that  (3.3)  is  well-posed  for  all  p.  assume  that  the 
input  z(x,p)  is  identically  zero  .  Now’  m(x,p)  is  the  linear 
least>squares  estimate  of  m(x,p)  based  on  observations 
which  are  identically  zero,  so  m(x,p)  =  Em{x,p)  =  0,  the 
last  equality  following  from  the  well-posedness  of  (3.1). 
As  a  result  X(x,p)  satisfies 

-£-\(x,p)  =  -A*\(x.p) 

and 

X(0.p)  =  \{L  j,p)  =  0 
and  thus  X(x.p)  =0  . 

4.  Recursive  Solution  of  the  Estimate  Hamiltonian 

We  shall  derive  our  recursive  smoothing  formula  by 
diagonalizing  the  dynamics  in  (3.3a)  via  a  change  of 
variables.  Let 

[  AP  <]BB' 

BP  =  r-\c>c  _  A‘ 

The  characteristic  equation  for  Hp  is  given  by 
X1  2Rc(o)X"  -f-  |a  |-  +  r~  ]q  =0 
wln're  a  —  p~-k2.  Solutions  to  this  equation  are 


The  four  eigenvalues  of  Hp  are  all  distinct  because  both 
the  real  and  imaginary  parts  of  (4.1)  are  non-zero  for  all 
p.  As  a  result,  we  can  diagonalize  Hp  as  follows  : 

,  PVr  0  1 

M~lHpMp  —  0  _Ap 

where 


Ap  [o  X0  ' 

and  the  modal  matrix  A/p  and  its  inverse  are  given  expli- 
citlv  bv 


Up  - 


X 0  X0  ~  ^0  ^0 

—  XqC  —  X0c  X q (i 


'<{  d\0~l  x0-‘  -r 

jq  -  c  -  c\0  1  -  X0  1 
•la  d  -  d\$  1  -  X0“  1  -  1 


-  r  c\r 


r  =  j(  ,)~o  )jq 


d  =  C  <7)/q 


a  —  (  qr~  1  2 


Now  l>v  using  the  following  change  of  variables  in  (3.3): 


losses*  aaaB»P5eaBB» 


,  /  I,P  =  Af-1  W’P) 

*»(^.P)  p  X(x,p) 


,'e  get 


d  fa  /  1  _  [A  p  0 
Tx  *b  =  0  -A 


VP  U  #y  Bf 
0  -A,  *,  +  Bb 


where 


Vf  Vj°]=  W0Mp  ,  \vf  Vjf-]  =  Wi.Mp 


(4.2) 


(4.3a) 


BI  =M-'  0 

B„  M'  -  C'r~ 1 


Equations  (4.3)  are  in  a  form  which  can  be  solved  recur¬ 
sively.  In  terms  of  'f/(0,p)and  4 'b(Lx,p),  a  solution  to 
(4.3a)  is 


where 


*/(s.p)  =  eXf{x)*  f(0,p)  +  9}(x,p) 

(4.4a) 

(x,p)  =  eA^Ll~z)^  b(Lvp)  +  $?(x,p) 

(4.4b) 

-~-V  °f{i ,p)  =  A p^J{x,p)  +Bfz(x,p) 

(4.5a) 

■^-4'A°(x,p)  =  -  Ap*?(x,p)  +Btz(x,p) 

(4.5b) 

*f(0,p)  =  0,  tf6°(Llfp)  =0. 

(4.5c) 

Note  that  (4.5a)  and  (4.5b)  are  stable  in  the  forward  and 
backward  directions,  respectively.  Setting  x  =  Lx  in 


,«V»V 


ft 


(4.4a),  and  x=  0  in  (4.4b),  we  can  solve  for  f(0,p)  and 
tyb(Lvp)  in  (4.3b)  as 

i.p)  +  yf*Ko.,)i 


Ffi  =  \vf  +  Vf  eA'il  |  Vt  +  ^°eA'L‘j 


where 


A  solution  to  (4.3)  is  therefore  given  by 


^/(x,p)  ^ 


0  -  e 


"*/(*. P) 

*b(x’P) 


\,(Ll-x)  p‘(^/0(ii.p)  +  n°*?(o,p)} 


(4.6) 


The  Fp,  matrix  is  invertible  because 


I  0 

Fp  =  fyA/,  q  gA>ii  , 

where  is  the  invertible  matrix  associated  with  (3.3) 
being  well-posed: 

Fh  =  W0  +  WLeH’L'  . 

The  behavior  of  the  algorithm  as  p-+  oo  needs  to  be 
examined  further.  We  will  show  that  the  determinant  of 
Fp  does  not  vanish  as  p —  oo.  As  p  gets  large,  one  can 
ignore  the  exponential  terms  in  and  write 

'0n+XoC  ^12  ^12 

-c  -d  0  0 

Fn  ~ 

$21  $21  ^22"^"XoC  6  *n'\-\(yd 


c 


where  8  tJ  is  the  ij-th  entry  of  n,,'1.  Then 
det Ffb  ~  -  (c-  d)2det(n„_  *)-  c'd2(\0-  \0)2+cd(c-  d)(\0-  \o)(0n+d22) 
As  p — *  oo,  \q—  Xq — *  0,  so  that 

lim  det F ,b  =  4( qr~  1+??2)9_2det(ni;_  1)y^0 

p  —00 

It  is  shown  in  Appendix  B  that  under  realistic  energy 
assumptions  on  the  observed  images,  $ j(x,p)  and 
'I >b(x,p)  decay  fast  enough  as  p-> oo  to  ensure  that 
m(x,y)  £  D(A)  and  X  £  D(A*). 

To  summarize,  the  solution  procedure  for  solving 
the  estimator  equations  (3.2)  is  (1)  transform  the  obser¬ 
vations  into  the  p-domain,  (2)  compute  ^f(x,p)  and 
'f'6°(x,p)  using  (4.5),  (l^find  V f(x,p)  and  <f>b(x,p)  from 
(4.6),  (4)  compute  m(x,p)  using  (4.2),  (5)  inverse 
transform  m(x,p)  to  get  m(x,y). 

5.  Other  Boundary  Conditions 

The  smoothing  problem  for  the  2-D  Helmholtz 
equation  with  more  general  boundary  conditions  can  be 
handled  in  a  manner  very  similar  to  that  just  described. 
In  addition,  a  wavenumber  that  varies  in  the  y  direction 
can  be  studied.  The  boundary  conditions  we  consider  in 
this  section  are  conservative  on  the  ^-boundaries,  and, 
in  general,  dissipative  on  the  x-boundaries.  To  derive 
the  estimator  for  such  problems,  we  start  off  with  the 
same  2-D  Helmholtz  equation  (2.1),  and  x-boundary 
conditions  (2.2a).  The  y-boundary  conditions  are  now 

cosij  u(x.O)+sin  d«(x.0)-0  (5-la) 
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cos-7  u(x,L  2)+sin7  uy{x,L2)  =  0  (5.1b) 

where  0  and  7  are  real.  Periodic  boundary  conditions 

u(x,0)  =  u(x,Lo)  ,  uy{x,0)  =  uy(x,L2)  (5-2) 

could  also  be  assumed.  The  wavenumber  k  satisfies 
k~{y)  =  k$(y)+jn,  r}> 0 

The  estimate  Hamiltonian  for  these  boundary  conditions 
and  wavenumber  is  the  same  as  (3.2),  however  the 
operator  T  is  defined  as 

V  =-{dyy+k\y))f 

with  domain 

D{  T)  =  {f  6  L 2(n  ):/ ,/y  abs.  cont ,  fyy  €  L2(fl  )  , 
f  satisfies  (5.1)  or  (5.2)} 

We  now  introduce  the  self  ad  joint  Sturm-Liouville  opera¬ 
tor 


Qf  =-(dyy+kg(y))f 


with  domain 


D(Q)  =D(T) 

The  operator  Q  has  a  countably  infinite  number  of 
eigenvalues  np  and  eigenfunctions  <J>p(y).  If  we  use  the 
transform  operator  K  defined  by 
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V  =  f*?(y)dy 


1 

i 


ihfcl  /W.  >  ^ 1  fcA 


From  this  point  on,  the  calculations  are  identical  to  the 
zero-boundary  condition  case,  except  that  K  replaces  S. 


To  consider  problems  where  the  x  boundary  condi¬ 
tions  are  more  general  than  (2.2a),  we  assume  that  the 
matrices  V0  and  VL  are  now  linear  operators  V0  and  VL 
such  that  under  the  transform  K  described  above,  the 
action  of  the  operators  V0  and  VL  are  multiplicative,  that 
is,  the  following  transform  relations  hold 

V0m(0,y)  <#=>  Vr0(p)m(0,p) 

VL  m(  L  v y )  <£=>  VL  ( p ) m(L  up ) 

where  F0(p)and VL(p)  are  complex  valued  2x2  matrix 
functions  of  p.  As  before,  it  is  necessary  that 
V0(p)  +  VL(p)eA,L  be  invertible  for  all  p.  An  example  of 
a  dissipative  boundary  condition  occurs  when 


s  1 

O 

lO 

V0(p)  = 

0  0 

.  VL(p)  = 

s  -  1 

and  s  is  complex.  This  case  corresponds  to  a  damped, 
elastically-braced  membrane  on  the  .x=0  and  x  =  Lx 
sides.  The  determinant  of  V0+  VL  eA,Ll  is 


Typically  (5.5)  is  non-zero. 


6.  Complexity  Analysis 

For  computational  purposes  the  rectangle  f]  is 
discretized  into  an  N  X  Af  grid.  The  complexity  of  the 
principal  steps  needed  to  calculate  the  leastrsquares  esti¬ 
mate  of  the  wave  amplitude  at  one  wavenumber  k  is 
given  in  Table  1. 


TABLE  1  The  Smoothing  Complexity  for  the  2-D  Helmholtz 
Equation 


Complex  it 

V 

Step 

zero  or  periodic 
v-bouncf  conds 

other  y- 

b.c. 

Find  (y)  and 
for 

(3.3U5.4) 

n/a 

0{M~b) 

Fourier 

transform  the 
observations 

0(NM\ogM) 

0{NM2) 

Find  m(x,p ) 

V  i<i 

(4.5), (4. 6). ( 4.2) 

O(NM) 

0(!VM) 

Inverse 

transform 

0(NM\ogM) 

0{NM2) 

In  Table  1,  b  is  the  bandwidth  of  the  matrix  used  to 
approximate  the  operator  Q  in  the  eigenvalue  calcula¬ 
tions.  Typically  b  —  l.  The  overall  complexity  for  each 
wavenumber  is  then  either  0(NM\ogM)  or  0(.VA/-). 
depending  on  the  y-boundary  conditions.  Obviously  if 
one  does  not  exploit  the  particular  structure  of  our 
model,  and  uses  only  first  and  second  moment  informa¬ 
tion,  the  complexity  would  be  0(M3N3).  As  discussed  in 
the  Introduction,  for  zero  boundary  conditions  a  com- 


plexity  of  0(.\/.Vlog.\/.V)  has  been  achieved  in  [2], [3], [4] 
using  different  techniques.  The  complexity  of  transform¬ 
ing  the  time  domain  images  into  the  frequency  domain 
for  T  wavenumbers  is  0(  TiXM log T).  Therefore,  the 
entire  smoothing  procedure  would  have  a  complexity  of 
0{TNM\ogTXf)  or  0{ TNMlogT+TNM2)  depending  on  the 
v-boundary  conditions  and  assuming  that  6=1. 


7.  The  Smoothing  Error  Covariance 

Using  results  in  [l],  we  can  express  the  smoothing 
error  rh(x,y)  —  m(x,y)~  m(x,y)  as  the  solution  of  the 
Hamiltonian  system 

£_\n(x.y)  1  [  -4  1BB']\m(x.y)  ]  \B  0  ]-  ,  x 

Ox  ~\{x,y)  r  lC'C  -A*  -\(x.y)  0  C'r  1  [u- (  x ,  y )  ] 


m  €  D(A )  ,  X  e  D(A  *) 

with  boundary  conditions  (compare  with  (3.2)) 

'  r  -  -i  p 

Un-'dv)  =  r'n - 1  v  +  0  /  7 

"  [y}  v  m(  L  ,.v)  0  /  -  \(  L  ,.v)  1  j 


Transforming  (7.1)  with  respect  to  y  gives  (see  Section 


' 1  m(  x  .p) 
:>x  X  ( x .  p ) 


Ir'CC 


qBU1  f  -  ,  .1  \b  0  lr  ,  \  n 

m(x,p)  u  c(J-p) 

-A’  -X(x.p)  0  C'r1  [u>(x.p) 


i-*n  -  i  t  \  _  it-  m(0,p)  ...  ™(L\’P)  t  -  oh 

l  n„  v(p)  o  -  X(0,p )  V/'  -X(Lj.p) 


Our  original  statistical  assumptions  and  the  properties  of 


iWy!*!1!1 


transform 


the 


imply 


th  at 


t(x,p)  ,  t(x.q)  ,  w{i,p)  ,  w{x,q)  ,  v(p)  ,  v(q)  are  mutually 
uncorrelated  for  pj^q.  Therefore,  if  we  let 


P(x,y)  =  £[m(i,y)m'(x,y)] 
P{x,p)  =  £’[m(x,p)m*(x,p)] 

then 


P(x,y)  =  E  P( x  <p)sm~py  ,  p  = 

Using  a  Green’s  function  solution  to  the  acausal  linear 
system  (7.2).  one  can  show  that  P{x.p)  satisfies 


where 


P(x,p)  =  [/  ;0]A(x,p) 


P) 


t, 

f  G{x, a  ,p) 
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0 

r-‘C'C 


G  *(x,<7  ,p)da 


+  ef{'IFH-'Vnv-lV(FH'y 


(7.3) 


and 


G(x,<7,p) 


X  >  <7 

x  <  a 


The  overall  complexity  of  solving  these  equations  is 
0(NM2),  due  to  the  integration  necessary  to  compute 
A (x,p)  in  (7.3).  An  alternate  procedure  with  complexity 
0( NM)  can  be  derived  by  first  diagonalizing  (7.2)  in  a 
manner  identical  to  that  in  Section  3.4.  If  we  change 
variables  using 


ef(x'Pi  ,f-  i  rh{  x  ,p) 
eb(x,p)  *  v  -\(x ,p) 


l  J  L  J 

then 


a 

e/{i,p) 

Ap 

0  ' 

e/(*,P) 

+- 

*/' 

p(*.p)  ‘ 

dx 

(X’P) 

0 

“V 

e6(*.p) 

5/ 

Mx'P) . 

i"n -lv(P) 
where 


v?  n° 


e/(0.p) 

e4  (O.p  ) 


tf{LuP) 
cb(L  i>p) 


15/ 

5/ 


M  - 1 
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0 
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A  solution  to  these  equations  is 
.  1  r  A ,  z  q 

'pj  =  C0  a,^,-*)  fy»l{^n;Wp)-^«/(^i.p)-  W(o,p)} 


.  /W) 

[e6°(^,P) 


where 


-ke°U'p)  =  +  B/  [*(  *’,'*)] 

(7.5a) 

-ie,0(x,p)  =  -A  ,e,V,p)  +  B!  [«/( x'.  pj  ] 

(7.5b) 

e/°(0,p)  =  0  ,  eb°{Lvp )  =  0 

(7.5c) 

We  can  now  express  A(x,p)  as 

A(x,p)  =  MpG{x,p)Mp  * 

(7.6) 

where 


©(j.p)  =  E  Wj’jjje/V-P)  eA 1  ■  P ) ) 


s 


We  will  express  the  diagonalized  error  covariance  ©(x,p) 
in  terms  of  the  second  moments  of  the  random  variables 
{v,ef(x),eb°(x)},  where  we  have  suppressed  the  argument 
p  for  simplicity.  Define  the  covariances 

Rf(xuXo)  —  E[e°(x1)ef°'(x2)] 

R/>(x  i,x2)  =  8le°(x1)eb°'(x2)} 

R/t(x  lfx2)  =  E[ef(x  i)ej0‘(ar2)] 

Now  0(x,p )  can  be  written  as 

*fl(z)Ffl  'fcf-  v'n;  'v+vfn,(i1)vj--+  vfR^L^o)  vf+ 


n°n»(o)  vr+vfRfrL  i,o)  vfX^)-'^(x) 

P/(x)  0 

-*fl(*)FfllG(x)-G ’(*)(Fjl)-l*jl(x)+  „  j,i(t) 

where 


eA,z  0 


*fb(X)=  q  eA  ,(£,-*) 


G(x)  =  Vf[Rf(Lltx)  !^y,(Ii,x)]fV*0[^(x,0)  ii?t(0,x)] 
Rf{Lx,x)  =  eA'(£'1‘z)n/(x)  ,  Rb(0,x)  =  eA,rn$(x) 
#/&(xl>x2)  =  -  eA'(l,~z*}nyj(x,-X2)  ,  X[>X2 


I/(x)  =Apn;(x)  +n;(x)A /+ 2170/0  °r  s/*’  n/(°)=° 


i»(x)  =Apnjx)  +n,(x)A/  +  °W*,  n,(L,)  =  o 


«.aTW 


n/6(X)  =  ^i7J0e  0  rJ^V'“< 


8.  Concluding  Remarks 

In  this  paper,  we  have  constructed  a  recursive 
smoother  for  the  2-D  Helmholtz  equation.  The  advan¬ 
tage  of  our  method  is  that  it  leads  to  computationally 
attractive  algorithms  even  when  the  x-boundary  condi¬ 
tions  are  dissipative  and  include  random  inputs. 
Another  advantage  is  that  a  linear  least^squares  Born 
inversion  of  the  wavefield  can  be  performed  along  with 
the  image  restoration  without  changing  the  computa¬ 
tional  complexity.  Higher  order  equations  characterizing 
vibrating  plates,  etc.  could  be  approached  in  a  manner 
similar  to  that  presented  here.  If  the  wavenumber  k  has 
variations  in  the  x-direction,  then  an  approach  based  on 
operator  Riccati  equations  could  be  used  to  formally 
diagonalize  the  smoother  dynamics.  However,  the  initial 
value  problem 

-7—  m{x.y)  =  Am{x.y) 
ox 

m{0.y)  =  m0(  ij) 

is  not  well-posed  fit  does  not  lead  to  a  semigroup),  and 
this  raises  questions  about  the  existence  and  uniqueness 
of  solutions  to  the  corresponding  Riccati  equations. 
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Appendix  A 

In  this  appendix,  we  show  that  the  state  variable 
representation  (3.1)  of  the  Helmholtz  equation  is  well- 
posed.  To  do  this,  we  show  that  the  transformed  system 
(3.4)  is  well-posed,  and  that  the  solutions  to  these 
acausal  linear  systems  give  rise  to  Fourier  coefficients 
that  ensure  that  the  formal  Fourier  series  does  converge 
and  that  the  state  vector  m(x.y)  e  D{A). 

(3.4)  are  well  posed  for  every  p  because  the  matrix 


V0  +  eA’L'  = 


1  0 
coshJL[  3~  ls\nh{3L  i) 


( A.l 


where  3  =  (p2-  k'2)1  is  invertible  for  all  p.  If  we  now 
change  variables  in  (3.4)  as  follows: 


T f(x,p) 
T  »(*./>) 


=  V 


m^x.p) 
m.i(  x  ,p  ) 


(A. 2) 


then 


Dr 


1  1 

-3  3 


IM 


d  T f(x.p) 
dx  T«fx .p) 


3  0 
0  3 


T f(x.p) 

T  b(x.p) 


+  1/2 


-3l 

3~' 


(^•P) 


A  .3  a) 
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-  21  - 


1  1  ‘ 

T  /(O.pf 
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o  o' 

T/fiLP)' 

0  0 

T6(0,p) 

i 

i  i 

T,(i,.p) 

v(p)  = 

We  can  write  the  solution  to  (A. 3)  as 


(A.  3b) 


T/(x.p) 

=  (l-e-2^1)"1 

0 

1 

_e-,Lr 

«i(p)-  TftO.p) 

0  'z) 

-e-^> 

1 

c2(p)-T«(Lltp)_ 

+ 


Tf(x,p) 

T4°(i,p) 


( A.l) 


where 


— T f{x,p)  =  -0Tf(x,p)  -  1/2  0  xe(x,p) 


-jjT 6°(x,p)  =  /?T4°(x,p )  +  1/2  0  le(x ,p) 


rf(0,p)  =T  ?(LllP)  =0 

Tf(x,p)  and  T6 (x,p)  are  continuous  functions  of  x  that 
depend  continuously  on  v  and  e.  Equations  (3.1)  will 
have  a  unique  solution  which  depends  continuously  on  v 
and  e  if  the  appropriate  Fourier  series  converge  and  if 
m(x,y)  e  D(A).  We  will  assume  that 

{d2/d'2y)v(y)  e  L22(0,Z/1)  ,  e(x,y)  6  L2(fl  )  and  show  that 
this  implies  that  ml(x,y)eD(T)  and  m2(x,y)  €  L2(fi ). 
We  first  derive  bounds  for  T °(x,p)  and  T6°(x,p).  Solving 
for  T f(x.p)  gives 


T?( x,p)  =  -±r f  e'**-'h(u,p)du 

-j  o 


which  gives 
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J  |e-  ^(jr_  “)  prfw 

/  N«,p)  l2rf« 
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/ S'  >SSSS 


'-Kv 
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/  \e-Kz-y>\2dy  < 


2  Re{0) 


Thus 


|T?(x,p)  |  < 


•11/2 


2\0  \V2Re(0) 

Since  0  «  p  when  p-> oo 

|T f(x,p)\2  ~  0(  Gp~3) 

where 


/  |e(u,p)  \2du 
o 


L  1 

G  =  f  \e(u,p)  1 2du 
o 

In  a  similar  manner, 

|T,°(i,p)P  -  0(Gp~ 3) 

Our  assumption  on  the  input  field  e  implies  that 
I/2 

/  ?dy<  00 

o 

almost  everywhere  with  respect  to  x  (a.e.  wrt  x  ).  By 
Parse val’s  theorem 


J]|e(x,p)  |2  <  oo  a.e.  wrt  x 
p 

and  hence  by  the  comparison  test  for  series 

|e(x,p)  |2  ~  o(p_1) 
a.e.  wrt  x  .  Therefore 


|T/(x,p)  |2  -  |Ta°(x,p)  j2  ~  0(p  4> 

The  first  term  in  (A. 4)  involving  the  boundary  term  v(p) 

is  becoming  exponentially  small  as  p-*  oo  for  all  x  in 

( 0, Z/ j ) .  Since  (d2/d2y)v(y)  e  L22(0 ,LX)  by  assumption,  the 

first  term  in  (A. 4)  has  a  convergent  Fourier  series  that  is 


twice  differentiable  with  respect  to  y  on  [O.Z.J,  and  is 
infinitely  differentiable  in  (O.Lx).  The  second  term  of 
(A. 4)  therefore  dominates  the  solution  in  (o.L,)  as  p 
gets  large,  hence 

|T/(*,p)|2  -  |T4(*,p)!2  -  0(p-<) 
a.e.  wrt  x.  By  using  the  change  of  variables  (A.l),  we 

see  that 

Imffx.p)  |2  =  \rf{x,p)+rb(x,p)  |2  -  o(p*4) 

and 

|m2(z,p)  |2  =  |-/3T/(x,p)-h/?T4(x,p)  |2  ~  0(p-2) 

These  bounds  show  that 

m  i(x,p)  G  D(T )  and  m2(x,y)  G  Lo(Q  )  [5] . 
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Appendix  B 

We  show  in  this  appendix  that  for  the  2-D 
Helmholtz  equation,  when  the  x-boundary  conditions 
are  separable  (0l2  =  92l  =  o  ),  and  the  observed  image  is 
bounded  and  mean-square  differentiable  in  the  in¬ 
direction, 

rii(x.y)  e  D{A)  and  \(x,y)  G  D{A‘) 

To  begin  with,  it  is  sufficient  to  show  that 
(<92/d2y)4'y (x,p)  and  (d2/d2p)'l»  4  (x.y)  exist  in  the  mean- 
square  sense.  The  boundary  conditions  will  be  satisfied 
due  to  the  sine  transform. 

The  p-domain  equations  for  V  f  4 , °. vp  b2°  are 


mmm 


given  by  (^/.-(x.y)  denotes  the  ith  component  of 
tyf{x,y),  etc.) 

X 

=  /  eHl-a)j—l - z(s,p)ds 

0  4A0r<7 

^/o°(x,p  =  -  f  =£ - z(s,p)ds 

0  4A0rcr 

Lx 

'i'  6  r°( x  ,p)  =  /  t~x^x~3)j~ q- - z{s,p)ds 

z  AXqTCT 

L  i  _ 

*bo°{x,p)  =  -  j  e~X^z~s)j—J- - z(s,p)ds 

z  4 \0ra 

As  in  Appendix  A,  one  can  show  that 

<  9ip~3/2M(p)  (B.l) 

where 

Lx 

M2{p)  =  /  | z{x,p)  1 2dx 
o 

and  gv  is  a  finite  constant.  The  same  result  holds  for  'ij. 
For  separable  boundary  conditions  [l] 

V  f{x,p)  =  *  f{x,p)+V  J{x,p) 

Where  for  large  p 

*fB(x,p)  =  [eA^  )0  - 1 K60^  6°(0.p  ) 

The  L2(Q  )  norm  of  V f(x,y)  can  be  expressed  using 
Parse val’s  relation  as 

oo 

l*/(x.ir)lfMn)  =  S 

p  =  1  0 


In  a  similar  fashion: 
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d2 


lx 


ll^*/(^)«V(n)  =  E  P4  J\\*f{x,p)\l2R*dx  (B.2) 

a“y  p  =  10 

We  wish  to  prove  that  this  norm  is  finite.  For  separable 
boundary  conditions  one  can  verify  that 

Vy'vft  < 


and 


|  |eA ' x  ;0 


<  9 3e 


ReX0i 


Combining  these  bounds  gives 
so  that 


lh^*/S(*,s)rLW  <  E  P3J,*W2(p)/e2R'x"'dx 

o  y  p  =  i  o 

or 

<  g  p-g^MHp) 

since  X0  <  g6p.  If  we  assume  that  ( d/dy)z(x,y )  e  L2(fi  ) 
then 


f)  00 

E  r«2(p)  o o 

y  p  =  i 

which  shows  that  ^ f[ x,y )  is  twice  differentiable  in  the 
y-direction.  Similarly, 

$2  oo 

l4-*°(*.p)Mn)  <  E  LipMHp)  <  oo 
which  implies  that  V f(x,y)  has  a  mean-square  second 
derivative  in  the  y-direction.  A  similar  argument  shows 
that  'f/b(x,y)  is  also  mean-square  twice  differentiable  in 
the  y-direction. 
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